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Abstract. We have investigated the quantum J\ — J2 — J3 model on the honeycomb lattice with exact 
diagonalizations and linear spin- wave calculations for selected values of J2/J1, J3/J1 and antiferromagnetic 
(Ji > 0) or ferromagnetic (Ji < 0) nearest neighbor interactions. We found a variety of quantum effects: 
"order by disorder" selection of a Neel ordered ground-state, good candidates for non-classical ground- 
states with dimer long range order or spin-liquid like. The purely antiferromagnetic Heisenberg model is 
confirmed to be Neel ordered. Comparing these results with those observed on the square and triangular 
lattices, we enumerate some conjectures on the nature of the quantum phases in the isotropic models. 
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1 Introduction 

Frustrated quantum antiferromagnetic (AF) spin systems 
on low dimension (D) lattices have attracted a great deal 
of interest in recent years. Quantum fluctuations, largest 
for small values of the spin S of the magnetic ions, low 
D and small coordination of the lattice, are expected to 
lead to novel magnetic behaviors. Their effects, have been 
preeminently seen in ID. They have been investigated on 
a few 2D systems. The most studied models are the AF 
Heisenberg model on the triangular or kagome lat- 
tice j3j which are geometrically frustrated systems, the 
J\ — J2 model on the square lattice QJ^,^| where frustra- 
tion is introduced by 2nd neighbor interaction, the Ji — J? 
model |?],|||| and the multi-spin exchange model (MSE) jio] 
on the triangular lattice. 

Less studied |ll|[l^ , p^| , ^4[ , spin models on the honey- 
comb lattice deserve attention due to the special proper- 
, ties of the lattice and since there are experimental realiza- 
tions. A first feature of the lattice is that, like the square 
lattice, it is not geometrically frustrated for AF nearest 
neighbor interactions but has lower coordination. Thus 
quantum fluctuations are expected to be larger than for 
the square lattice. For this reason the spin- 1/2 Heisenberg 
antiferromagnet on the honeycomb lattice, has been stud- 
ied theoretically by various methods Jllp^jr^ , p^ | which 
all predicted that Neel long range order (LRO jsubsists 
but with an order parameter smaller than for the square 
lattice case. This also motivated a Schwinger-boson study 
of the effects of frustrating second neighbor interactions 
in the J\ — J2 model [jl5| . 

A major incentive to study frustrated magnets on the 
honeycomb lattice is the availability of experimental data 



in the family of compounds BaM 2 (X0 4 ) 2 (M= Co, Ni; X= 
P, As) obtained, some years ago, by Regnault and Rossat- 
Mignod jlfj]. The magnetic ions M have small spins (it 
is supposed to be S = 1/2 for the Co oxide and 5 = 1 
for Ni), disposed in weakly coupled layers where they sit 
on a honeycomb lattice. The simplest model relevant to 
these quasi 2D compounds is a J\ — Ji — J3 model on a 
honeycomb lattice with first, second and third neighbor 
interactions and either on site if (S = 1) or XXZ (if S = 
1/2) anisotropy. 

So far the J\ — J2 — J3 model was only investigated 
within first order linear spin-wave theory (LSW) 0,0, 
and to our knowledge the renormalization of the order 
parameter by quantum fluctuations has not be calculated 
even in this simplest approach. The experimental results 
motivated us to do this calculation in the large S limit and 
then attack the 5=1/2 problem with exact diagonaliza- 
tions (ED). 

In this paper, as a first step, we consider the case 
of purely isotropic interactions. The Hamiltonian of the 
model reads: 

H = Ji Si-Sj+Ja S - S ^ + J 3 Y S - Sfe W 



<i,j>l 



<i,k> 2 



<i,k> 3 



where the first, second and third sums run on the first, 
second and third neighbor pairs of spins, respectively (see 
Fig.|l|). The coupling constants Ji can be either AF (Jj > 
0) or ferromagnetic (Ji < 0). Depending of the values 
of the parameters Ji, this model displays various classi- 
cal ground-states: a collinear AF ground-state, two de- 
generate manifolds of planar hclimagnctic ground-states 
with four or eight sublattices and a ferromagnetic ground- 
state. The classical phase diagrams of the isotropic and 
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Jj>0 




J 2 A 



Fig. 1. The honeycomb lattice. The full and empty circles 
differentiate the two sublattices. ti and t2 are the two vectors 
of the triangular Bravais lattice. Top right: arrows show the 
interactions between a site and its 1th, 2th and 3th neighbors. 



anisotropic XXZ models display only minor differences. 
In particular the classical ground-states are the same for 
the parameters believed to be relevant to BaCo2(As04)2- 
In this paper we concentrate on the quantum effects in the 
isotropic model. The study of the quantum XXZ J\ — J2 — 
J3 model with parameters appropriate to BaCo2(As04)2 
will be presented in a separate paper fl8|| . 

We investigated quantum effects for selected values of 
the Ji chosen to give a broad picture of the quantum ef- 
fects encountered in the model. The restriction to isotropic 
interactions limits the number of parameters and enable 
to separate the effects of anisotropy. In addition to the 
J\ — J2 models on the square and triangular lattices, the 
present model may be compared with the J\ — J2 — J3 
model on the square lattice which has a similar variety 
of classical ground-states and to which a few studies have 
been devoted |l|£||P||]2|§| . 

This paper is divided into five parts. In section II, we 
recall the classical phase diagram of the model, obtained 
by Rastelli el al. |17| , and identify the degeneracies of the 
classical ground-states not considered by these authors, we 
also discuss the stability of the first order spin-wave ap- 
proximation for these different phases. The ED results for 
the case of antiferromagnetic and ferromagnetic nearest 
neighbor coupling are presented in section III and IV re- 
spectively. In section V we draw conclusions, and enumer- 
ate some conjectures relative to the appearance of the var- 
ious generic two-dimensional spin-liquids. We described in 
an appendix the various technical features specific to our 
present ED calculations on different samples. 



Fig. 2. Classical phase diagram for antiferromagnetic near- 
est neighbor interactions. In the T = classical approximation 
regions II and IV have in fact a degenerate manifold of non- 
planar ground-states. Thermal fluctuations or quantum fluc- 
tuations do select the collinear configurations shown in this 
figure. 

2 Classical phase diagram and semi-classical 
deviations 

2.1 Planar ground-state configurations 

The classical model was studied by Rastelli et al. 
They searched for planar or uniformly canted configura- 
tions minimizing the classical energy E c i . The former were 
found energetically favored over the latter. They represent 
spiral configurations, characterized by a wave- vector Q. 
The classical spin (of length S) sitting at cell R of the 
triangular Bravais lattice on sublattice a is given by: 



Sr, q = S^cos (Q.R + <p a ) u + sin (Q.R + 



(2) 



where u and v are two orthogonal unit vectors defining 
the plane of the spins, 4> a can be chosen to be zero on one 
sublattice and will be noted <f> on the other. The set of 
spiral wave- vectors Q minimizing the classical energy will 
be noted {Q}. 

The phase diagram of planar solutions of this type is 
reproduced in Fig]| (Ji > 0) and Fig.|| (Ji < 0). There is 
a mapping between the two phase diagrams: the transfor- 
mation J\ — > — Ji, J3 — > — J3 and S,; —>■ — Si for i £ on the 
black triangular sublattice of Fig.l leaves the Hamiltonian 
unchanged, and maps the ground-state for J\ > on that 
for J x < 0. 

There are six regions in each phase diagram: four collinear 
phases (I, II, IV, VI) and two spiral ones (III,V). In the 
collinear regions I and VI, the wave- vector of the magnetic 
order is Q = 0, whereas in regions II and IV, {Q} = {Ki}, 
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Fig. 3. Classical phase diagram for ferromagnetic nearest 
neighbor interactions. Same comments on regions II and IV as 
in Fig.| 



where are the three inequivalent middles of edges of 
the Brillouin zone (see Fig]|). The phases are: In I, <f> = it 
(0) if Ji > ( Ji < 0), in VI, <j) = (tt) if Ji > (Ji < 0), 
in II <j> = ir (0) if Ji > (Ji < 0), in IV 4> = (tt) if 
Ji > (Ji < 0). 

The separation lines I-III, I-V, II-III, IV-V represent 
continuous transitions; the others are first order phase 
transitions. 



2.2 Non planar ground-states manifolds 

Ansatz (Eq.||) is usually generic to find all allowed classi- 
cal ground-states |26|j27[ . It assumes that, up to the triv- 
ial degeneracy associated to a global spin rotation, the 
ground-state is unique or exhibits at most a discrete de- 
generacy. It is valid if a linear combination of the different 
Q modes of the same set {Q}, can be excluded as violat- 
ing the constraint |Sr jQ ,| = S on every site. Exceptions 
occur for special sets {Q}, in particular if Q is half or 
one fourth of a reciprocal lattice vector G [^6| . This is the 
case in region II and IV where {Q} = {G/2} = {Ki} (see 
Fig. 6). Here, the linear combination of three Kj solutions 



>R., 



2_, S cos (Kj.R + <p a ) Ui 



(3) 



i=l 



with unnormalized u^, is submitted to the constraints u^.Uj 
Si j and uf = 1. The ground-state is a two dimensional 
manifold continuously connecting the three solutions 
(there are nine degrees of freedom for choosing the three 
Uj, minus three for global rotations of the spins and four 
constraints). This gives birth to the non planar ground- 
states manifolds described below. 



Fig. 4. Top: four-sublattice classical ground-state in region IV 
on Fig.^ for antiferromagnetic first neighbor coupling (Ji > 0). 
Bottom: the collinear solutions with the three possible arrange- 
ments (in this case, classical spins in sublattices A and B are 
antiparallel). 



In regions IV for Ji > or II for J\ < ( (j> = 0), the 
ground-state manifold is the set of four-sublattice ordered 
solutions such as Sa + Sb + Sc + Sd = 0. This could be 
seen directly from the expression of the classical energy of 
these configurations: 



Ed = jj(Ji + 2J 2 ) {Sa + S b + S c 
- ^ ( Ji + 2 J 2 - 3 J 3 ) (S* + S| 



(4) 



si,) 



In this equation, N represents the total number of 
spins of the sample. The generic four-sublattice config- 
urations are shown in Fig^f, as well as the three collinear 
configurations, which appear as special cases of it, with 
Sa = Sb = — Sc = — Sd and the two other combi- 
nations of parallel spins. The situation is reminiscent of 
the Ji — J2 model on the triangular lattice with 
1/8 < J 2 /Ji < 1. 

In regions IV for J\ < or in II for Ji > (<fr — 
7r) , the ground-state manifold is the set of eight-sublattice 
solutions shown in Fig.|| where sublattices labelled by the 
same letter are paired (partner sublattice is over-headed 
by a bar) and 

Sg = S a (5) 
for a G {A, B, C, D} and 



S A + S B + S C + S D = 0. 
This minimizes the classical energy: 



(6) 



N 



J 2 (S, 
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Fig. 5. Same as Fig.^) but for ferromagnetic first neighbor 
coupling (Ji < 0), region IV of Fig|j. 

+S^ + S B + S + S D ) 2 

+ ^(Ji- 2^2) (3 A + S b + S c + S d ) 

(S^ + + S<5 + S 5 ) 
+ ^ (3J 3 - Ji) [(S A + S A f + (S B + S s f 

+ (S C + Scf + (S D + S B ) 2 } 
- 1 (- Ji + 2J 2 + 3 J 3 ) (Si + S| + S^ + S£ 

+S| + S| + S^ + S 2 5 ). (7) 

It is highly probable that thermal fluctuations will sta- 
bilize the collinear solutions as it does in similar situations 
on the square and triangular lattice. We will show below 
that quantum fluctuations indeed do it. 

Continuous degeneracy of the ground-state also occurs 
when Q = G/4 in III and in V but since this happens 
only on lines (for instance if J 2 = 0.5 in III) and not 
in full regions, we shall skip their description which is 
not essential to our present goal. The transition line III-V 
between the two spiral regions is very special. It has an 
infinite degeneracy of spiral ground-states corresponding 
to: 

cos(Q.ti) + cos(Q.t 2 ) + cos(Q.(ti - t 2 )) = 1/8 Jf - 3/2 

(8) 

The lines of Q solutions of Eq.|| are shown in Fig.|| for 
J 2 = 0.2,0.4,0.5. 

2.3 Stability of the quasi-classical phase diagram in 
the large S limit: LSW results 

The renormalization of the order parameter m) in the first 
order spin-wave approximation is already large in the pure 



Fig. 6. Brillouin zone of the triangular Bravais lattice. Light 
solid, dashed and dotted lines are the solutions of Eq.8, for 
J2 = 0.2, 0.4, 0.5 respectively. 



model (J 1 — 1, J 2 = 0, J3 = 0): w) w 0.48 i.e. a value 
reduced to w 48% of its classical value (FigjTj). 

For Ji > 0, the interplay of quantum fluctuations and 
frustration quickly destroys Neel LRO: m) goes to zero 
for J 2 w 0.1 or J3 w —0.1. The helimagnetic LRO of zone 
V disappears: near the J3 = axis, this is mainly due 
to the large classical degeneracy of the ground-state and 
near the J 2 = axis, the main cause is the vanishing of the 
spin-wave velocity at the point J\ = 1, J 2 = 0, J3 = 0.25. 
The zone V being very small, we conclude that the Neel 
helimagnetic ground-state does not survive in region V for 
antifcrromagnetic J\. 

For Ji < the ground-state is ferromagnetic in zone I 
of figure The ferromagnetic state is an exact eigenstate 
of the hamiltonian, quantum fluctuations don't destroy 
it. However, the classical degeneracy on the boundary be- 
tween zones III and V implies a whole branch of soft modes 
and a disappearance of the helimagnetic LRO, on and near 
this line (as in the AF J\ case) . The main difference with 
the AF first neighbor case is the persistence of Neel order 
in zone V near the J\ = —1, J 2 = 0, J3 = 0.25 point and 
in the vicinity of the J 2 = axis. 

The nature of the quantum spin- 1/2 phase for small J 2 
and J3 appears an open problem that we will now attack 
with the help of exact diagonalizations (ED). 

Exact diagonalizations were performed on samples of 
N = 18,24,26,28,30,32 sites with appropriate bound- 
ary conditions (see Appendix). The technical problems 
encountered in such approaches have already been stud- 
ied in details in previous references ||§ and will not be 
described in details in this paper. We briefly discuss in Ap- 
pendix the different characteristics of the studied samples 
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Fig. 7. LSW values (triangles) of the order parameter rrv for 
the purely AF Heisenberg model (Ji = 1, J2 = 0,J3 = 0), as 
a function of N 1 ^ 2 . The asymptotic behavior ~ N 1 ^ 2 (dashed 
line fitted to the LSW values) is only reached for quite large 
samples, much larger than the sizes studied in exact diagonal- 
izations. 

with respect to the present model. We will now proceed 
to the analysis of the ED results. 

3 J\ > 0: AF nearest neighbor interactions 

3.1 The purely Heisenberg model, and phase I of the 
quantum model 

The AF Heisenberg point {J\ = 1, J% = J3 = 0) was previ- 
ously investigated by exact diagonalization calculations on 
small samples |ll]] , Monte Carlo simulations |Q , series ex- 
pansions around the Ising limit p3| , spin- wave theory up 
to second order |l4| , and Schwinger-boson mean field the- 
ory ]Ti|. All concluded that the quantum system exhibits 
Neel LRO but with a large reduction of the order parame- 
ter due to quantum fluctuations (the Monte Carlo result, 
rather close to the spin- wave value, is m) = 0.44 ± 0.06). 
Our ED results for sizes up to N = 32 are consistent with 
this conclusion. The approach is however different. In an- 
tiferromagnets with linear Goldstone modes, the scaling 
law for ru* is an expansion in 1/N 1 / 2 [ [29|]3Cfl . For the 
sizes encountered in ED, the asymptotic 1/N 1 / 2 law is 
never reached and the extrapolation of the order pa- 
rameter to the thermodynamic limit remains uncertain. 
A qualitative idea of the finite size scaling of rrv can be 
obtained from the LSW results: they show that the asymp- 
totic l/N 1 ! 2 regime cannot be expected for sizes smaller 
than N~ 400 (see Fig. 0). 

Nevertheless confirmation of Neel LRO can be obtained 
thanks to characteristic features of the spectrum itself 
which have more favorable scaling behaviors: 



Fig. 8. AF Heisenberg model, scaling of the QDJS with S 
and N for N = 18, 24, 26, 28, 32. 



For a given sample size, the lowest eigenlevels in each 
sector of the total spin S evolve as Eq(S, N) oc S(S + 
1) up to 5 ~ V~N, as shown in Fig|| They are the 
eigenlevels associated with the collective dynamics of 
the order parameter, the so-called Quasi Degenerate 
Joint States (QDJS)[|lJ, which can be described by the 
effective Hamiltonian: 



Hcoll.dyn. — —A(N) S' 



(9) 



where A{N) is the finite-size difference in total energy 
between the absolute ground-state and the first triplet 
excitation. 

The number of QDJS and their symmetries are those 
expected for the projections of the classical Neel order 
on the irreducible representations (IR) of 5(7(2) CS> Q 
({?:lattice symmetry group): There is just one state 
for each S value since there is just one way to couple 
two spins of magnitude A74 in a total spin S. These 
states are invariant under lattice translations (their 
wave- vector is k = 0), under 2-7r/3 rotations around the 
center of an hexagon, they are even (odd) under inver- 
sion with respect to the center of an hexagon for even 
S and N — 4p (respectively odd S and N — 4p + 2). 
In the thermodynamic limit, the QDJS collapse on the 
singlet ground-state as 1/N (Fig. |l0|b). The QDJS re- 
main distinct from the softest magnon excitation which 
collapse to the ground-state as 1/N 1 / 2 . The QDJS and 
the softest magnon are shown for the = 32 sample 
in Fig|. 

The asymptotic ~ 1/iV behavior of A(N) is not yet 
reached for our largest sizes as seen in Fig.^ojb. The 
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Fig. 9. Low energy part of the AF Heisenberg spectrum for 
N = 32: eigenenergies are plotted versus eigenvalues of S 2 . Full 
triangles represent QDJS; empty squares describe the softest 
magnon. The dashed-line and the dotted-line are guides to the 
eye for the QDJS and the softest magnon respectively. 



Fig. 10. AF Heisenberg model, (a) energy per site eo: the line 
is a fit to the leading term of eq.ll. (b) spin-gap: the full line 
is a fit to eq.10, the dashed line if a linear fit in 1/N. 



3.2 Region IV 



next order term in the 1/N 1 / 2 expansion reads f2{||3(| : 

4(,v) ^ (1 -"^ ) + 0( ^ ) (10) 

where \ is the spin susceptibility, c is the spin-wave 
velocity, p the spin stiffness and (3 is a number of order 
one. A fit of the spin-gaps to this scaling law is shown 
in Figjl(i|b. The importance of the term in 1/N 3 1 2 is 
not unexpected since this term is oc c/p and quantum 
fluctuations which reduce p with respect to its classical 
value are strong (as already shown by the reduction of 
the order parameter). 
— The energy per site eo(N) — Eq(0, N)/N of the ground- 
state scales as: 

for a Neel order. FigjlOja shows that the leading term of 
order 0(1/N 3 ^ 2 ) is enough to describe the size effects 
in this range. 

A rapid analysis of the quantum phase diagram in re- 
gion I does not reveal new phases, but both LSW cal- 
culations and ED confirm that a weak antiferromagnetic 
second or ferromagnetic third neighbor coupling are suf- 
ficient to kill LRO: the boundary between phase I and 
phases III and V is shifted upwards by quantum effects. 



In region IV the classical model presents a degenerate 
manifold of four-sublattice ordered ground-states. The finite- 
size spectra clearly show that this degeneracy is lifted by 
quantum fluctuations which favor a collinear two-sublattice 
order (see Fig.[ll]): the low lying levels of these spectra, be- 
low the magnons excitations, exhibit a large family { 4 E} 
of QDJS associated to four-sublattice solutions. At the 
bottom of this family there appears a line of cigcnlcvels 
with definite symmetries: these levels constitute the family 
{ 2 E} of QDJS states associated to a collinear symmetry 
breaking. 

The situation, seen here, is very similar to the one 
previously studied for the J± — J2 antiferromagnet on the 
triangular lattice ||. The expected number of states 4 Ns, 
2 Ns in { 4 E} and { 2 E} and their space symmetries are 
easily determined for each value of the total spin S. The 
eigenstates of { 4 E} can be labelled by the five irreducible 
representations i~j (i — 1, 5) of S4 (permutation group of 
four elements) . The mapping between the space group op- 
erations on the four-sublattice solutions and permutations 
of 54 is described in Table E], together with the charac- 
ter Table of S4. The four-sublattice order is invariant in 
two- fold rotations (TZ^): thus the eigenstates of { 4 E} be- 
long to the trivial representation of C2 . Since it is also in- 
variant under a two-step translation of the Bravais lattice 
they have either a wave- vector k = or a wave- vector . 
Ji , J2 and J3 belong to the k = subspace, whereas P4 
and I5 belong to the subspace {K^}. r± and J2 are invari- 
ant under the three- fold rotations 7?. 27 r/3 of C3, whereas 
P3 is associated with the two-dimensional representation 
of C% v . The number of replicas of Ji that should appear for 
each S value can be computed as in ref || . The result are 
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Fig. 11. Low energy spectrum for Ji = 1, J2 = 0, J3 = —2 
and iV = 32. Full triangles represent states of the family { 2 E}, 
empty triangles represent states belonging to { 4 E} and not 
to { 2 E}, and empty squares represent the softest magnon. All 
theses states have the symmetries predicted in Tables 2 and 3. 



Fig. 12. Ji = 1, Ji = 0, J3 = —2, energy gaps measured 
from the absolute ground-state versus 1/N for N = 24, 32. Full 
squares connected by the dashed line: gap to the lowest energy 
state in the triplet sector (it belongs to { 2 E}). Full triangles: 
gap to the 2 nd singlet state of symmetry i~3 (it belongs to 
{ 2 E}). Open triangles: gap to the 3 rd singlet state of symmetry 
i~3 (this state belongs to { 4 E} and not to { 2 E}). 



shown in Table g for the N=32 sample. Analysis of the 
two-sublattice order can be done similarly: the collincar 
solution has a three-fold degeneracy, the set of eigenstates 
{ 2 E} maps on Z 3 . It is characterized by the IR TV -T 3 and 
7^4 . The number of repliquas are shown in Table || for the 
N=32 sample. 

The "order out of disorder" phenomenon pq| is clearly 
seen for J\ = 1, J 2 = 0, J 3 = —2. In Fig]ll| we show 
the lower part of the N — 32 spectrum at this point. 
The lowest eigenstates in each S sector are the states of 
{ 2 E}, describing collinear order. Further support to this 
assumption is given by the finite size effects of the energy 
gaps. As shown in Fig.[l2| a plot of the spin-gap of the N = 
24,32 samples versus 1/N is consistent with a vanishing 
value for N — * 00. On the other hand the gap between the 
two states i~i and I3 of { 2 E} of the 5 = sector tends 
to close when the size goes to infinity, whereas the gaps 
between the levels of { 2 E} and the other levels of { 4 E} 
increase with N. 

We have investigated the scaling behavior of the spec- 
tra at some other points of region IV not too close from 
the classical boundaries and found essentially the same 
behavior and a selection of collinear LRO by quantum 
fluctuations. 

Closer to the boundary between region IV and V, the 
separation between the { 4 E} states and the magnons states 
decreases. This is an indication of a softening of the magnons 
and the neighborhood of a 2nd order phase transition to- 
wards another phase. The behavior of the spin-gaps at 
Ji = 1, J 2 = 0JL J 3 = -0.5 and J x = 1, J 2 = 0, J 3 = -1, 
similar to FigJlJ, nevertheless indicates that these points 



of the quantum phase diagram are still in the collinear 
phase IV. 

Various studies of the spectra of the N = 18, 26, 28 
samples under suitable boundary conditions confirm these 
results for the quantum phase IV, and indicate that the 
quantum boundary between phase IV and V is probably 
slightly shifted down relatively to the classical boundary 
shown in Fig. 2. 

Results of ED calculations (not shown) in region II for 
Ji < (which has the same classical manifold of degen- 
erate ground-states as in IV for J x > 0) suggest a similar 
selection of the collinear solution there too. In conclusion 
up to a slight motion of the boundaries , the semi classical 
behavior in regions I, II, IV and VI, is not qualitatively 
affected by the strong quantum fluctuations of the spins 
1/2. 



3.3 Quantum phases between I and IV 



The intermediate phases between the two collinear Neel 
phases cover region V and part of region III. In this part 
of the quantum phase diagram, SU (2) symmetry is unbro- 
ken, and there is a gap to triplet excitations: these phases 
only support short range order in the spin-spin correla- 
tions (see Figjl3|). Our LSW and ED calculations indicate 
that this quantum region likely extends in regions I and 
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Fig. 13. Spin-spin correlations as a function of distance for the 
pure Heisenberg model on N=24 (triangles) and N=32 samples 
(squares) and for Ji = 1, J 2 = 0.3 on N=24 (three-legged star) 
and N=32 samples (four-legged star). 



IV [] In this work we study region V, region III close 
to I, and the transition line III-V with ED calculations 
using TBC on N = 18, 24 samples (see Appendix) and 
PBC on the N = 24 and N = 32 samples. A thorough 
search of the ED spectra, sweeping the twist angles at the 
sample boundaries, did not yield evidence of incommen- 
surate helical LRO, neither with the wave- vectors of the 
classical solutions nor at other wave-vectors. In all cases 
no tower of QDJS was found. The ED results corrobo- 
rate the conclusion of LSW calculations that the classical 
spiral solutions are destabilized by quantum fluctuations. 
This seems a rather general statement in systems where 
the quantum fluctuations are strong enough [po|j2l]| . 

Is this quantum phase a quantum disordered one? To 
answer this question we performed extended ED calcu- 
lations on N = 24, 32 samples on different points of the 
transition line III-V where the classical model has an infi- 
nite set of spiral ground-states, and the LSW calculation 
diverges. Along this line, we found evidence of two differ- 
ent phases both with a gap. 

Let us begin by the phase around J2 = 0.4: this point 
is very close to the point where the energy versus J2 is 
the largest and may be considered as a point of maximum 
frustration. The spectrum of the N = 32 sample is shown 
in Fig.[l4| This spectrum differs from the spectra of the 
collinear ordered system in IV: 



the lowest states are not IR of { 2 E} 



1 LSW calculations predict non vanishing order parameters 
for the classical spiral solutions inside III away from the bound- 
aries but a vanishing order parameter in the whole region V. 



Fig. 14. Low energy spectrum for Ji = 1, J 2 = 0.4, J3 = and 
N — 32. Full triangle: ground state; empty triangle: first singlet 
excited state (these two states have a wave vector k = 0); 
empty square: second singlet excited state; full square: first 
triplet state. 



and features associated with Neel LRO are missing: 

— The lowest eigen-energies for each S value do not in- 
crease as S(S + 1) with S 

— The lowest states in each S sector are not separated 
from the others as the QDJS are separated from the 
magnons, instead there is a dense continuum of states 
in each S sector except the S — one. 

— Furthermore a plot of the spin-gap of the N = 24, 32 
samples versus 1/N, displayed in Fig.[l5|, shows that 
the scaling law characteristic of a Neel ordered system 
is not obeyed and indicates a large spin-gap ps 0.2 for 
N -> 00. 

Most likely however the system is not fully disordered 
but exhibits dimer LRO (see FigJT^). The dimer opera- 
tor on a pair of sites is d{j = (1 — Pij) /2 where 
Pjj = 2(S,.Sj + 1/4) is the spin permutation operator. 
This projector is greater (less) than 0.25 when the spin- 
spin correlation is negative (positive), equal to 1 on a 
singlet and to on a triplet. For J\ = 1, J2 = 0.4, 
on the N = 32 sample, the first neighbor correlation 
is < dk,i >= 0.4899. On the symmetry breaking Spin- 
Peierls state (pure product of ordered dimers), the aver- 
age value of the dimer operator is 1 on the bonds where 
there is a dimer, and 1/4 on the other bonds. As the 
exact eigenstate does not break C3 symmetry the num- 
ber 0.4899 should be compared with the result obtained 
on the symmetric superposition of the three Spin-Peierls 
states aligned along the three main directions of the lat- 
tice: in this symmetrized Spin-Pcicrls state this correlation 



is d 



0.5. The average value of the dimer operator 
in the exact ground-state is thus very close to the Spin- 
Peierls value. 
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Fig. 15. Ji = 1, J2 = 0.4, energy gaps measured from the 
absolute ground-state versus 1/iV for N = 24, 32. Full squares: 
spin gap, i.e. gap to the first triplet excitation; open triangles: 
gap to the first singlet excitation (k = 0, IR F^); open squares: 
gap to the second singlet excitation. 



The dimer-dimer correlation between a reference bond 
and the bond (k, I) is — < di,jdk,i > — < 

dij >< d)-,i >• As in ref fiofl , we normalized -D(i,jV(fc,r) 
by its maximum value which is achieved when the two 
bonds are completely correlated. We thus measured dimer 
correlations by 



D 



P(i,j),(k,l) 



< dk.i > - < dij X dk,i > 

< dj,jd k ,i > - < djj >< d k ,i > 

< d hj > (1- < d k: i >) 



(12) 



If pujY(k,l) = 1 the presence of a dimer on bond im- 
plies the existence of a dimer on bond (A;, I); if Pu,j),(k,l) 
there is an absence of correlations between singlets on 
bonds and (k,l). If Pu,j) (k,t] ^ s negative, a singlet 

on bond induces a tendency towards ferromagnetic 
correlation on bond (k, I). 

The correlation pattern for dimer on first neighbor 
bonds is displayed in Fig.|lq. The calcul of dimer-dimer 
correlations on the state ^f y ^ gives P(ij).(k,i) = +0.5 if 
(ij) and (k,l) are parallel and Pu,j),(k,l) = —0.25 if (i,j) 
and (k,l) are non parallel bonds. The exact dimer-dimer 
correlations are not too far from these values and decay 
very slowly with distance. This is in favor of a columnar 
LRO of dimers with a C3 symmetry breaking, previously 
proposed by Einarsson and coll. Q . 

Moreover the degeneracy of the ground-state for N — ► 
00 points to the same conclusion: Figjl^ indicates that the 
gap between the and I3, S = 0) lowest states, which 
are both k = states, closes for N — » 00, while the 
gap between these states and the upper levels increases 



Fig. 16. Ji = 1, J2 = 0.4, singlet-singlet correlations 
P(i,2),(fe,i) between the reference bond (1,2) and bonds (k,l) 
in the ground-state of the N = 32 sample. The numbers above 
bonds (k,l) are the values of P(i 2),(*,i) truncated to the two 
first significant digits. Full (dashed) lines indicate positive (neg- 
ative) values of P(i,2),(fc,o , The width of the lines is proportional 
to the magnitude of |p(i,2),(k,i) I- The number above the bond 
(1,2) is < di,2 > (see text). 



with the size. Ji is non degenerate and J3 twice degener- 
ate: this allows the building of the three columnar dimer 
patterns with a C3 symmetry breaking and no transla- 
tion breaking. In this picture the finite size ground-state 
is the symmetric combination of these three states. This 
degeneracy corresponds to a true symmetry breaking with 
a local non zero order parameter (dimer LRO): this is a 
Valence-Bond Crystal(VBC). 

The honeycomb lattice could a priori accommodate a 
different kind of VBC with alternation of hexagons with 
three dimers and hexagons without dimers: this pattern 
breaks both C3 and translational symmetry. In their large 
N approach, Read and Sachdev |Uj found that this struc- 
ture might be the ground-state. In the range Ji « 0.3 — 
0.35, we find a short range structure roughly reminis- 
cent of this arrangement. In fact the short range dimer- 
dimer correlations are even more symmetric than in this 
VBC crystal and would be more compatible with a crys- 
tal of hexagon-plaquettes in a symmetric 5 = state |p3| : 
for example the correlation (1-2) (7-6) should be negative 
and equal to —0.25 in the Read and Sachdev VBC state 
whereas it is equal to 0.1 in the pure hexagon-plaquette 
VBC. In fact in our SU(2) model all correlations decrease 
noticeably with distance (Fig. [l3[ ^) and the pattern does 
not seem to propagate at large distances. The ground state 
in this range of parameter is probably a RVB spin liq- 
uid. This conclusion is qualitatively substantiated by the 
study of the energy gaps to the ground-state: plausibly 
none of them goes to zero at the thermodynamic limit, 
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Fig. 17. Singlet-singlet correlations for Ji = 1, J2 = 0.3 and 
N = 24 (same legend as in Fig.|l^) . 



which would be consistent with the RVB hypothesis. Un- 
fortunately the finite-size effects are rather chaotic: the 
ground-state energy of the N = 18 and N = 30 samples 
(samples on which the Read and Sachdev dimer pattern is 
not frustrated) are larger than the ground-state energy of 
the N — 26 and N — 32 samples, which frustrate it. This 
is probably related to the fact that the N — 18 and N = 30 
samples do not allow the system to take full advantage of 
the second neighbor antiferromagnetic coupling, whereas 
the N = 26 and N = 32 samples do. But the building of 
singlets on second neighbor bonds tends to destroy VBC 
patterns and favor a RVB ground-state. All these argu- 
ments point in favor of an RVB phase at this coupling 
parameter: unfortunately the sizes that can be studied do 
not allow a quantitative determination of the gaps. 

The quantum AF J\ — J2 model on the square lat- 
tice close to the point of maximum frustration exhibits 
aimer § or plaquette f§ LRO; the same kind of conclusion 
has been drawn for the J1 — J3 model S| , and for the MSE 
model on the square lattice |25j . These phases share qual- 
itative properties with the phase identified for J2 = 0.4: 
in each cases a collinear LRO is destabilized by frustra- 
tion giving birth through a 2nd-ordcr phase transition to a 
massive phase with dimer LRO. Such VBC phases appear 
in many models on bipartite lattices: Rokhsar and Kivel- 
son |34j in the Quantum Dimer approach (QD), Dom- 
bre and Kotliar [|35) for the Hubbard model, Read and 
Sachdev [^J in the SU(N) approach of the Heisenberg 
model found VBC phases. These phases, as the first one 
described here (for J 2 = 0.4), have a gap for all excita- 
tions, a discrete degeneracy of the ground-state, exponen- 
tial decrease of the two points spin-spin correlations but 
LRO in higher correlation functions; they have only con- 
fined spinons. 



Fig. 18. Low energy spectrum for Ji = — 1, J2 = 0.25, J3 = 
and N = 32. Full triangle: ground state; empty triangle: first 
singlet excited state (these two states have a wave vector k = 
0); empty square: second singlet excited state; full square: first 
triplet state. 



Up to now we only know few spin- 1/2 models exhibit- 
ing true RVB phases with a clear-cut gap: the MSE model 



on the triangular lattice [ 1 
on the triangular lattice] 



and the Quantum Dimer model 
J . More work is needed to know 
if the excitations of these different RVB states are simi- 
lar and in particular if they sustain deconfined spinons 
excitations. 



4 Ji < 0: ferromagnetic nearest neighbor 
interactions 



As already underlined above in Sect. 3.2 the classical collinear 



phases (F or AF) observed for large Ji and J 3 are likely to 
survive to quantum fluctuations. We thus focus our study 
on the region of maximum frustration, < | J%\, \ J$\ < 0.5, 
corresponding to region V and part of region III of FigJ^. 
In this situation LSW calculations predict a non vanish- 
ing order parameter of the spiral solutions for values of 
J2 and J3 not too close to the transition lines. However 
extensive ED calculations performed on the N = 18, 24 
samples with twisted boundary conditions do not yield 
any evidence of spiral LRO. 

We thus studied samples up to N = 32 spins to in- 
vestigate the nature of the ground-state for a few sets of 
parameters. The most extensive calculations were done at 
the point J2 = .25, J3 = on the transition line III-V 
which may be considered as highly frustrated (for this J2 
value the ground-state energy is close to its maximum, 
and in the LSW approach quantum fluctuations destroy 
LRO). Strong indications that the model has a RVB spin- 
liquid ground-sate, were found at this point: 
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Fig. 19. Ji = — 1, J2 = 0.25, energy gaps measured from the 
absolute ground-state vs 1/N for TV = 24,32. Black squares 
show the spin-gap. Open triangles (squares) the gap to the 
first (second) excitation in the singlet sector. 



— The spectra do not exhibit a tower of QDJS as shown 
in Fig.|l8| for N = 32, and Eq(S) clearly does not evolve 
as S{S + 1) with S. 

— A plot of the spin-gap versus 1/N, shown in Fig.[l9|, 
indicates that the spin-gap is small but finite when 
TV — > 00 |. 

But contrary to the case with positive J\ and J2 = 0.4: 

— The correlations display a strong short range order but 
plausibly no LRO. The short range pattern is original: 
the first neighbor spin-spin correlation is ferromagnetic 
(< Si.Sj >= 0.10), the second (third) neighbor spin- 
spin correlations are antiferromagnetic ( < Si.Sj > n . n .= 
—0.13, < Si.Sj > n .n.n.= — .25 ), but no long range 
pattern does emerge from this picture. The dimer- 
dimer correlations equally show a strong short range 
pattern and apparently no LRO. Fig. g(J represents 

2 In view of Figjl9| one may object to our extrapolation 
to TV — » 00 on two numerical samples. In fact our conclu- 
sion is supported by examination of both the gap and the 
energy per bond of the samples with sizes 18, 24, 28, 32 
with various boundary conditions (available on request at: 
fouet@lptl.jussieu.fr). This study shows that the variations of 
these quantities with the size is very small and mainly due 
to the frustration of the short range antiferromagnetic order 
between third neighbors due the boundary conditions (see be- 
low) and not to the cut-off in the low-energy long wave-length 
quantum fluctuations. Notice that the energy per bond on the 
two non frustrating sizes 24 and 32 does not increase with the 
system size but decreases by a very small amount (~ 10 -3 ). 
We thus conclude that we are in the cross-over regime for both 
sizes 24 and 32 and the extrapolation of the spin gap in Fig.|l9| 
is reasonable. 



Fig. 20. Ji = — 1, J2 = 0.25, first neighbor triplet-triplet 
correlations on the N = 32 sample. The spin-spin correlation 
on the reference bond (1, 2) is ferromagnetic: the number above 
this bond measures the projection of the two-spin state of the 
exact ground-state on the pure triplet state. 

first-neighbor dimer-dimer correlations: they are much 
weaker than in the AF case (remark that triplet-triplet 
correlations are equal to sing let-singlet ones and com- 
pare with Fig. [if]). Fig. £l] displays second- neighbor 
dimer-dimer correlations which also decrease with dis- 
tance. The third-neighbor dimer-dimer correlations de- 
crease even quicker. The strength of the short range 
correlations explains the finite size results on small 
frustrating sizes (see footnote 4 ). 
— The ground-state is probably unique in the thermody- 
namic limit. The first two singlet excitations are shown 
together with the first triplet excitation in Fig. |l9|. The 
spin-gap and the ground-state energy per spin display 
a very small sensitivity to the size for N = 24, 32. The 
gap to the third excitation seems more sensitive to the 
size but this might due to the fact that the different 
sizes do not accommodate the same wave vectors. In 
view of the results it seems probable that the singlet 
excitations will not collapse to the absolute ground- 
state in the thermodynamic limit. 

We thus conjecture that the quantum ground-state of 
this system does not break any symmetries of the Hamilto- 
nian or of the lattice: it is a quantum spin-liquid where all 
excitations are gapped. Results obtained for the N = 24 
and 32 samples for J2 = .5 also indicate a finite spin-gap 
when N — > 00. This suggest that there is a spin- liquid 
phase in a finite range of parameters. 

Such a quantum massive phase, without LRO, is highly 
reminiscent of the spin-liquid phase found in the MSE 
model on the triangular lattice |n^. Curiously enough it 
appears in the two cases in the vicinity of a ferromagnetic 
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Fig. 21. Ji = — 1, J2 = 0.25, second neighbor singlet- singlet 
correlations on the N = 32 sample. The reference bond is 
(1, 10): the spin-spin correlation on this bond is antiferromag- 
netic, the number above the bond measures the projection of 
the two spin-state of the exact ground-state on a pure singlet 
i.e.: < di.10 >; the value of this observable in a symmetrized 
wave function of products of second-neighbor singlets is 0.375. 
Non-parallel singlet-singlet correlations have been omitted for 
clarity, all of them are very small and negative. 

phase destabilized by antiferromagnetic couplings. There 
is a difference in the degeneracy of the ground-state in the 
two cases: whereas the ground-state on the honeycomb lat- 
tice is unique, it has a 4-fold degeneracy on the triangular 
lattice. This is easily understood as the honeycomb lat- 
tice is not a Bravais lattice and has two spin- 1/2 per unit 
cell. Thus the uniqueness of the ground-state in this latter 
case does not contradict the Lieb-Schultz-Matthis- Affleck 
conjecture pq |, or the topological approach of Read and 
ChakrabortypT] . 



5 Conclusions and conjectures 

This study of the spin- 1/2 J\ — J2 — J3 model on the 
honeycomb lattice has brought the following new results: 

— For small frustrations J2/J1 or J3/J1 less than ~ 0.15 
or larger than 1, the system remains essentially clas- 
sical: when various kinds of LRO are possible, quan- 
tum fluctuations, as well as thermal fluctuations in the 
classical case, select the LRO with the most symmetric 
order parameter amongst the various possibilities. 

— The classical symmetry between the phase diagram for 
ferromagnetic J\ and the phase diag ram for antifer- 
romagnetic J 1 discussed in Sect. |2.l| is destroyed by 
quantum fluctuations. 

— For the largest frustrations these models exhibit gapped 
phases. 

— For an antiferromagnetic first neighbor coupling, a Va- 
lence Bond Crystal phase has been clearly evidenced 
around J 2 / J\ — 0.4. 



— For an intermediate frustration J2/J1 — 0.3, an RVB 
spin-liquid appears between the Neel ordered phase 
and the VBC phase. 

— For a ferromagnetic first neighbor coupling, the present 
results favor the hypothesis of a RVB spin-liquid phase 
in a large range of parameters. No VBC has been found 
in that case. 

This study of the spin-1/2 J\ — J2 — J3 model on the 
honeycomb lattice, when compared to similar approaches 
of SU(2) Hamiltonians leads us to formulate some conjec- 
tures on the generic behavior of such models on different 
lattices. 

— In 2D the pure S=l/2 Heisenberg model is Neel or- 
dered on any bipartite lattices with coordination num- 
ber > 3. It is disordered on the triangular-based kagome 
lattice which has a coordination number equal to 4. 

— Non-coplanar classical ground-states are not robust 
against quantum fluctuations in the isotropic models. 

— Neel order or ferromagnetism disappears around the 
points of maximum classical frustration giving birth 
to phases with spin-gap and short range spin-spin cor- 
relations. 

— Disappearance of a ferromagnetic phase due to anti- 
ferromagnetic frustrations leads generically to a spin- 
liquid phase, with short range correlations in all ob- 
servables and a gap to all excitations. 

— Disappearance of a collinear antiferromagnetic phase 
might lead to a VBC phase either directly (Ji — J2 
model on the square lattice), or through an interme- 
diate RVB phase (this study). The spin-liquid phase 
observed by Santoro et al in the spin-orbital model ]3S|] 
might be rather similar to the RVB phase described 
here. 

For completion we might add: 

— Disappearance of a non-collinear phase (3-sublatticc 
Neel phase) takes place through a phase with a spin 
gap but a continuum in the sing let sector 0,0. 
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Appendix: Special properties of the studied 
samples, and boundary conditions 

The ED calculations were performed as in refs |l],|l(| on 
systems of N = 18, 24, 26, 28, 30, 32 sites shown in Figg 
The N = 18, 24, 32 samples have the full point group sym- 
metry of the lattice, whereas the N = 26 sample misses 
axial but still has rotational C3 symmetry, the N = 28, 30 
have neither. 

With periodic boundary conditions (PBC), all the sam- 
ples are of course compatible with the Q = order in 
region I. In region IV, however only the N = 24 and 32 
samples have the full symmetry of the classical order. The 
N — 28 sample is compatible with one collinear solution 
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Table 1. Character table of the permutation group S4. First 
line indicates classes of permutations. Second line gives an el- 
ement of the space symmetry class corresponding to the class 
of permutation. These space symmetries are: t the one step 
translation (A — > B), 7^2vr/3 (^2^/3) the three- fold rotation 
around a site of the D (B)-sublattice, and a the axial symme- 
try keeping invariant C and D. N e i is the number of elements 
in each class. 



Fig. 22. The N = 18, 24, 26, 28, 30, 32 samples. 



but frustrates the two others as well as the non coplanar 
solutions. The other samples are frustrating but can allow 
a collinear order if twisted boundary conditions (TBC) are 
used. This is the case for the TV = 18 sample if a twist of 
7r is applied along ti or t2- 

To search for spiral order, we used TBC and sweep the 
whole range [0, 2ir] of twist angles (pi.2 in the ti and t2 
directions. These specific boundary conditions are defined 
as: 

S(R i +t j ) = e ^^( R -)s(R 4 )e-*^ (R ' ) . (13) 

This allows to look for boundary conditions which would 
not frustrate helical ground-states. This approach was found 
effective for the Heisenberg model on the triangular lat- 
tice to deal with samples frustrating the three-sublattice 
Neel order Q. For such samples the ground-state energy 
was found to reach its minimum for the twists which re- 
lease the frustration: at that point the spectra recover the 
characteristic features of Neel order. 

A VBC with the pattern of Read and Sachdev 
(considered in Sect. fits in the N — 30 sample but 
not on a N = 32 sample. 
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